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Abstract. We implement the damage spreading technique on 2-dimensional isotropic and anisotropic Bak- 
Sneppen models. Our extensive numerical simulations show that there exists a power-law sensitivity to the 
initial conditions at the statistically stationary state (self-organized critical state). Corresponding growth 
exponent a for the Hamming distance and the dynamical exponent z are calculated. These values allow 
us to observe a clear data collapse of the finite size scaling for both versions of the Bak-Sneppen model. 
Moreover, it is shown that the growth exponent of the distance in the isotropic and anisotropic Bak-Sneppen 
models is strongly afi^ected by the choice of the transient time. 



PACS. 05.65.-|-b Self-organized systems - 64.60.Ht 
evolution 

1 Introduction 

In 1993, Bak and Sneppen (BS) introduced a simple model 
to describe the biological evolution of an ecology of in- 
teracting species [1]. Since then, it has attracted quite 
some attention among the natural scientists as well as 
economists. This wide range of interest comes from the 
fact that BS model is the simplest model that exhibits 
self-organized criticality [T]f5]. The self-organized critical- 
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ynamic critical phenomena - 87. 23. Kg Dynamics of 

ity feature of the BS model is revealed in its ability to nat- 
urally evolve towards a scale invariance stationary state 
[3]. That is, the correlation length in the BS model is infi- 
nite and an initial local perturbation might lead to a global 
effect. Therefore, it is worth to study the sensitivity to the 
initial conditions in the BS model. 

The technique that we use to study the sensitivity to 
the initial conditions in the BS model is known as the 
damage spreading technique in dynamical systems theory 
and can be described as follows: If we consider two copies 
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of the same dynamical system starting from slightly dif- havior of the temporal and spatial correlation functions 

ferent initial conditions and follow their time evolution, with Id isotropic and anisotropic BS models. Furthermore, 

the sensitivity function can be defined as, we expect topologically the same short-time evolution of 

Ax(t) the Hamming distance with different scaling exponents. 

at)= Ijm — AZ^cxp(At), (1) 

In order to identify the dynamical quantities in 2d BS 

where A is the Lyapunov exponent and Ax(0) and Ax(t) , , -n ■ ; i • ,i ^ ,• ,i n 

model we will mtroduce m the second section the well 

are the distances between two copies at t = and i, re- , . . , . . . „ „„ , , 

fcnown isotropic and anisotropic versions oi the ab model 

spectively. Depending on the A being positive, negative , i; r ,i •,• i , i r 

and compare our results tor the critical exponent values of 

or zero three different behavior can be distinguished from , , , . , , ^ . , , n i • 

temporal and spatial correlations to the known nndmgs. 

Eq. (fTl): (i) A > 0, the system is said to be strongly sen- t .i ,i • i ■• • i ; ;i i t 

in the third section, we implement the damage spreading 

sitive to the initial conditions, (ii) A < 0, the system is , , . r. i • , • i • ; • -^.n ^ i i 

technique on Zd isotropic and anisotropic Bb models and 

said to be strongly insensitive to the initial conditions, , , , ,, i , i •,• , , mi i i 

calculate the related critical exponents, ihe relation be- 

(iii) A = 0, while the sensitivity function could be a whole , , , i i • • i 

tween some model parameters and the critical exponents 

class of functions, for the low-dimensional discrete dynam- . . , i- 

IS reviewed. Moreover, we analyze the hnite size scaling 

ical systems Eq. (Il]) forms a power-law, r .i i- i tt ■ i- ; <• i ;i • r 

J ~L r , normalized Hamming distance for both versions of 

£,{t) ~ t", (2) the BS model. Finally, in the last section, we discuss the 

results and give some perspectives. 

where a is some exponent and a > and a < cases 
correspond to weakly sensitive and weakly insensitive to 

the initial conditions, respectively [31[511M3- For the high 2 Isotropic and anisotropic versions of BS 

dimensional systems like the BS model, the same analysis model 
could be performed using the so-called Hamming distance 

instead of the sensitivity function. The dynamics of BS model on a 2-dimensional lattice of 

Our task will be to investigate the short-time dynam- an edge size L has simple rules. The initial state of the 

ics of the isotropic and anisotropic BS models on a square system on a square lattice is characterized hy N ^ L x L 

lattice using the standard damage spreading technique. As fitness values (random numbers) fij, where i = 1---L 

far as we know, this technique has not been used in the and j = I - ■ ■ L, uniformly distributed between and 1. 

literature to investigate the sensitivity to the initial con- These fitness values are assigned to each site i, j of the 

ditions in 2-dimensional {2d) BS model while there have lattice with periodic boundary conditions. The usual dy- 

been several works in Id [51l51 [TUl[TTlll2j . In 2d, with ex- namics of the system is eventually achieved by localizing 

tensive simulations, it is our hope to obtain similar be- the lattice site with the minimum fitness fmin and assign- 
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ing new random numbers to that site and its first nearest ulations since the convergence is faster than the interme- 

neighbors. diate anisotropy choices [13] . 

The model described above can be called as an isotropic , • ^ i- i - i j j • 

^ Alter some transient time which depends on the size 

BS model since the interaction between the current mini- r , i ^ xi ■ ^ • j • ^ • j i i • 

01 the system the isotropic and anisotropic models achieve 

mum and its first nearest neighbors is the same in both di- ^ ^- ^- ^^ ^ ^- ^ ^ i- ir • i 

^ a statistically stationary state (i.e., seli-organized criti- 

rections. In other words, if one observes the minimum fit- i ^ ^ \ • i • i i r ^ 

cal state) m which the density oi htness values is uni- 

ness value of the system at time step t as ^mi■a — fi i then r ij-4--u4-j r-rii j -t, m-fi u 

J f J nun iormly distributed on [jc, IJ and vanishes on [0, jcj, where 

the fitness values fi^ij, ftj, fi+i,^ and critical threshold value /c depends on the lattice di- 

will be updated at that time step. This means that the . ^ n i- • ^ ^ ■ ^ i i 

mension. Un the 2-dimensionai lattice considered here, 

possibility for the minimum to iump to one of its left or „ „ „„„ i„ ^ ■ ^ • no i i i j- « 

^ -It ^ 0.328 [3 for the isotropic Bb model and jc — 0.439 

right nearest neighbors or to one of its up or down nearest r ,i ■ , -no i i rm j-jx i , ,i 

° ° for the anisotropic BS model, ihe dmerence between the 

neighbors in the next time step is simply the same, with •,• i ,i i u i r • ^ • j • ^ • no 

° ± ± J 7 critical threshold values of isotropic and anisotropic Bb 

a resulting isotropic avalanche of events. , , r ,i i • ,i ^ c^r. j- 

° models comes from the change m the rates of the spreading 

On the other hand, one can easily consider alterna- out of avalanches. Once the stationary state is achieved, 
five updating rules. One possible alternative could be to the temporal and spatial correlation functions are power- 
update at each time step fi-u.j, fi.j-e, fi,j, fi+b,j and law in both models signifying the existence of a critical 
fi,j+r (where u, £, b and r are arbitrary positive integers state with no characteristic length or time scales (scale 
taken from the interval [1, L\). Let us consider, for exam- invariance). These correlations can be used to determine 
pie, the case of fi~2.j, fi,j-i, fi,j, /i+i, j a-nd /i, j+2. The the universality classes of such models. The distribution 
system now has an inherent bias to the up and right. This of the absolute distance x between successive minima is 
means that we would expect a preferred direction for an a good example for the spatial correlation and defined 
avalanche to propagate. Therefore, this model is called as by Pjump{x) ^ x^'^ , where tt ~ 2.92 for the isotropic 
the anisotropic BS model. In principle, several types of BS model and tt ~ 2.57 for the anisotropic BS model on 
anisotropy can be introduced by changing the values of u, a square lattice. The temporal correlations Pfirst{t), the 
i, b and r, provided that £ ^ r and u =/= b (since £ ~ r distribution of first return times, and Paii{t), the distribu- 
and u = b represents the isotropic BS model). The maxi- tions of all return times scale as Pfirst{t) t^'^ first g^j^^l 
mal anisotropic cases are defined by = (£ = 1), r = 1 Paii{t) ^ t^"^"" , where Tfirst — 1-24 and Taii — 0.71 for the 
(r = 0) and 6 = (6 = 1), u = 1 (u = 0) whereas other isotropic BS model; Tfirst — 1-32 and Taii ~ 0.85 for the 
definitions are considered as intermediate anisotropics. In anisotropic BS model. These different values of tt, Tfirst 
this work, we use the maximal anisotropic case in all sim- and Taii suggest that the isotropic and anisotropic BS 
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Fig. 1. Time evolution of the normalized Hamming distance for different transient times of the isotropic (left) and anisotropic 
(right) BS models on a square lattice with A*' = 200 x 200 fitness values. As the transient time increases, the slope value 
decreases and eventually, converges on a fixed value. For TV = 200 the slope values are a ~ 0.66 and ol ~ 0.74 for the isotropic 
and anisotropic BS models, respectively. 

models belong to different universality classes. Our simu- shown that this distance exhibits an initial power-law di- 
lation results for the critical exponent values of temporal vergence, followed by a finite size-dependent saturation 
and spatial correlations given above are in good agreement regime, 
with the known results for 2d BS model [51 [nill41ll5| . 



3 Damage spreading 

As mentioned in the Sec. IJ), the sensitivity to the initial 
conditions of the isotropic and anisotropic BS models on 
a square lattice can be studied by damage spreading tech- 
nique. This technique has been used in the literature pre- 
viously to investigate the propagation of local perturba- 
tions in \d BS model [SligifTUllTTIITClfTTlfT^ . In these works, 
measuring the evolution of the discrepancy between two 
initially close configurations under the same noise, it was 



The algorithm for the BS model on a square lattice 
with N = LxL fitness values can be introduced as follows: 

1. Once the stationary state has been achieved, consider 
the system as replica 1 denoted by f} ^ 

2. Produce an identical copy of f} ^ and introduce a small 
damage in this copy by interchanging the site with 
minimum fitness with a randomly chosen site (denote 
the new replica as ff^j)- 

3. Let both replicas [flj and ff j) evolve in time using 
always the same set of random numbers. 
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Fig. 2. Time evolution of the normalized Hamming distance for five different system sizes of the isotropic (left) and anisotropic 
(right) BS models. In ^ oo limit, the slope values are estimated as a = 0.83 ±0.03 for the isotropic (left) and a — 0.91 ±0.03 
for the anisotropic case (right), respectively (see Fig. ([3|). 



We can now define the Hamming distance between two 
replicas as, 



' -I 



(3) 



where (■ • ■ ) stays for the configurational averages over var- 
ious realizations. Although, one of us introduced a new 
definition of the Hamming distance for the damage spread- 
ing technique [TB1IT71ITS] , which allows one to analyze both 
short- and long-time dynamics, in this work we use the 
standard definition given in ([3]), since we aim to analyze 
the short-time dynamics better by performing extensive 
simulations. 

In our simulations, all the results arc obtained by av- 
eraging over 800 realizations, which is enough to avoid 
the fluctuations. The largest system size that we achieve 



is TV = 400 X 400 = 1.6 x 10^^ for both isotropic and 
anisotropic versions of the BS model. In order to let the 
system achieve the statistically stationary state we choose 
an appropriate transient time for each lattice size. The 
choice of the number of transients is important to obtain 
the most realistic slope value. For a fixed lattice size, one 
can attain a couple of different transient times each of 
which is enough to observe the power-law growth of the 
Hamming distance. But, as it can be seen from Fig. ([1]), 
the smallest transient time gives the largest exponent. As 
the transient time increases, the exponent value decreases 
until it reaches a fixed value and further increment of the 
transient time does not affect the slope anymore. This is 
the most reliable way of choosing the appropriate tran- 
sient time for each system size. Otherwise one can easily 
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Fig. 3. The confidence interval of a in ^ cxd limit. From this scaling we estimate the growth exponent of Hamming distance 
as a = 0.83 ± 0.03 and a — 0.91 ± 0.03 for the isotropic (■) and anisotropic (t) BS models, respectively. 

overestimate the exponent of the power-law growth. In our 
simulations, we choose this number as 4 x 10^, 7 x 10^, 
1 X 10^ 2.4 X 10^ and 6.4 x 10^ for 50 x 50, 75 x 75, 
100 X 100, 200 X 200 and 400 x 400 lattices, respectively. 
Once the stationary state is achieved by attaining the most 
appropriate transient time, one can study the short-time 
dynamics of the model using the measure given in Eq. ([3]) . 
For Id isotropic and anisotropic BS models, it is known 
from the literature that this measure exhibits initially a 
power-law growth such as -Do(0 ^ ^" saturates at a 



constant value [8l[9lfT0 j fTT | ll6 1 [T7 l fT2] . As it is evident from 
Fig. 12]), our simulation results claim that the procedure 
that describes the time evolution of the Hamming distance 
in Id isotropic and anisotropic BS models remains true 
also for 2d isotropic and anisotropic BS models. The lim- 
iting case exponents are obtained as a ~ 0.83 ± 0.03 and 
a ~ 0.91 ± 0.03 for the isotropic and anisotropic versions 
of the BS model, respectively using the extrapolation in 
Fig. ([3]). These results of 2c? can be compared to the ones 
coming from the Id cases (namely, a = 0.48 and a — 0.53 
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Fig. 4. Log-log plot of r versus for the five curves given in Fig. ([2]). The estimated dynamical exponent values are z = 0.84 
and z — 0.73 for the isotropic (■) and anisotropic (t) BS models, respectively. 

3.1 Dynamical exponent and finite size scaling 

The second important exponent defined by z is called the 
dynamical exponent and comes from the scaling of t{N) ^ 
, where r is defined to be the value of t at which the 
power-law increasing part of the Hamming distance mea- 
surement crosses over onto the saturation regime for fixed 
N (namely, in Fig. intersection of two straight lines 
drawn through the linearly increasing power-law curve and 
the horizontal constant plateau). From Fig. (jl]), we ob- 



for the Id isotropic and Id anisotropic cases, respectively 
|12j). The increase of the exponent as the dimension in- 
creases, for both isotropic and anisotropic cases, reveals 
that the ability of the system to cover the whole lattice 
increases with the dimension. It is also possible to compare 
these results with the mean-field values (~ 1) obtained in 
[10] for the ring model and in [19^ for the coherent noise 
model. 
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Fig. 5. Data collapse of finite size scaling given in Eq. (j4| for the five curves of Fig. ([2]) in the isotropic (left) and anisotropic 
(right) BS models. 



tained the dynamical exponent as 2 ~ 0.84 and z ~ 0.73 
for the isotropic and anisotropic versions, respectively. 

The obtained exponents a and z can be used to an- 
alyze the finite size scaling behavior of the isotropic and 
anisotropic BS models on a square lattice. Namely, us- 
ing the exponent values one can control if the normalized 
Hamming distance D{N,t) ~ {D{t) / D(l)) obeys the fol- 
lowing finite size scaling behavior. 



D{N,t)^NPF(-^ 



(4) 



where (3 — a"f and j ~ z {(3 ~ 0.69, 7 ~ 0.84 and (3 ~ 0.66, 
7 ~ 0.73 for the isotropic and anisotropic cases, respec- 
tively). Substitution of corresponding exponent values in 
Eq. ([U leads to a clear data collapse for each version of BS 
model, as it can be seen in Fig. ([5]) . Different values of the 
dynamical exponent obtained for isotropic and anisotropic 



cases as well as Id and 2d cases indicate that the dynam- 
ical property, characterized by the time needed to cover 
all lattice sites, differs for each case. These values are also 
different from the mean-field model results (~ 1) [TUIITO] . 
which clearly states that the time for the system to reach 
the plateau scales linearly with the lattice size. 

4 Discussion and conclusion 

We study the short-time dynamics of 2d isotropic and 
anisotropic BS models on a square lattice. Our exten- 
sive simulation results for the power-law growth of the 
temporal and spatial correlations are in good agreement 
with the known results of the previous works. Moreover, 
the time evolution of the Hamming distance is studied 
by implementing the damage spreading technique and it 
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is revealed that 2d isotropic and anisotropic BS models 
exhibit weak sensitivity to the initial conditions as it is 
the case for Id BS model on a chain. It is shown that 
the scaling of the Hamming distance strongly depends on 
the choice of the transient time. That is, if one chooses a 
relatively large transient time for a fixed system size N 
to have the system at the statistically stationary state, a 
power-law behavior of the Hamming distance in time can 
be observed. But it is seen that increasing the transient 
time leads to a smaller slope value a. The most realis- 
tic transient time can be attained by observing the slope 
value of the time evolution of the Hamming distance until 
it reaches a fixed value. Eventually, we find the values of 
the parameter a equal to 0.83 ± 0.03 and 0.91 ± 0.03 by 
implementing the damage spreading technique on the five 
different lattice sizes, 50 x 50, 75 x 75, 100 x 100, 200 x 200 
and 400 x 400, of the isotropic and anisotropic BS models. 
The numerically obtained a values can be considered as a 
good estimation for a in iV ^ oo limit and together with 
the dynamical exponent z lead to a data collapse of the 
finite size scaling for the five different lattice sizes of the 
isotropic and anisotropic BS models. 

The authors wish to thank L. F. Lcrnrncns for discussions 
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(Turkish Agency) under the Research Project number 104T148. 
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